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We present the results of the first study of global oscillations of relativistic stars with both 
elastic crusts and interpenetrating superfluid components. For simplicity, we focus on the 
axial quasi-normal modes. Our results demonstrate that the torsional crust modes are essen- 
tially unaffected by the coupling to the gravitational field. This is as expected since these 
oscillations are known to be weak gravitational- wave sources. In contrast, the presence of a 
loosely coupled superfluid neutron component in the crust can have a significant effect on the 
oscillation spectrum. We show that the entrainment between the superfluid and the crust 
nuclei is a key parameter in the problem. Our analysis highlights the need for a more detailed 
understanding of the coupled crust-superfluid at the microphysical level. Our numerical re- 
sults have, even though we have not considered magnetised stars, some relevance for efforts 
to carry out seismology based on quasi-periodic oscillations observed in the tails of magnetar 
flares. In particular, we argue that the sensitive dependence on the entrainment may have to 
be accounted for in attempts to match theoretical models to observational data. 
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1 Introduction 

Following the observations of quasi-periodic oscillations (QPOs) in the tails of giant flares in several soft 
gamma-ray repeaters (SGRs) [TJ EJ [3] , axial (torsional) oscillations of neutrons stars have attracted consid- 
erable interest. The SGRs are generally thought to be highly magnetised neutron stars (magnetars) [4], and 
it would be natural to explain QPOs in the range ~ 30 — 150 Hz as axial crustal oscillations. That magnetar 
activity might trigger such oscillations was, in fact, suggested quite some time ago [5]. This explanation is 
supported by the spacing of the QPO frequencies [5]. The model does, however, face a serious challenge 
since the strong magnetic field should couple the crust to the core in less than an oscillation period [5] . 
In addition, the magnetic field will affect the motion of the crust itself [9l QUI Qj]. To determine global 
"elasto-magnetic" oscillation modes is challenging due to the possible existence of an Alfven continuum in 
the core. The presence of this continuum casts doubt on the very existence of global mode solutions with a 
discrete frequency spectrum, see [121 113] for recent progress on the purely fluid problem. However, recently 
there has been some evidence |14( 115] in favour of the hypothesis put forward in [5] . For moderate magnetic 
fields (B < 10 15 G) the global modes most easily excited by a catastrophic event in the crust have frequencies 
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that are tuned to the purely elastic crustal mode frequencies. For very high magnetic fields ~ fO 16 G the 
situation is not quite as clear [IB] . 

Most available studies of the axial mode problem have completely ignored the dynamical role of the 
dripped neutrons in the inner crust. Since these neutrons are likely to be superfiuid they can flow through 
the crust lattice. This leads to the problem having an additional degree of freedom, and it is important to 
establish to what extent the presence of this new component will affect the seismology. A recent local analysis 
[TTj has shown that the superfiuid components, both in the crust and the core, may have significant effects. 
The main aim of this paper is to consider this problem in more detail. We will, for the first time, examine 
the effects of the superfiuid crust component on the global axial oscillation spectrum. The obtained results 
have immediate relevance for the magnetar discussion. We are also reporting real progress in modelling 
the dynamics of realistic neutron stars. Our discussion provides insight into the global dynamics of solids 
coexisting with a superfiuid component, and highlights the microphysics parameters that are needed if we 
want to improve the analysis. As we will see, the entrainment between the superfiuid neutrons and the crust 
nuclei (the "effective mass" of the free neutrons) is a key parameter that needs to be constrained better by 
equation of state calculations. 

Our analysis is based on state-of-the-art matter modelling. We employ a linear perturbation scheme in 
general relativity, and do not make additional approximations such as the Cowling approximation, wherein 
the gravitational degrees of freedom are neglected (or, as discussed in [BJ, the coupling is ignored). Since 
we account for the dynamics of spacetime itself, the oscillation modes we determine can be split into two 
distinct families, the torsional i-modes associated with the matter motion supported by the elastic properties 
of the solid, and the gravitational w-modes (which in turn can be subdivided into different classes, see e.g. 
[15]). In principle, we should be able to address directly the gravitational- wave damping of the t- modes and 
the influence of the elastic and superfiuid nature of the matter on the w-modes. In practice, however, we 
are unable to compute the damping time of the i-modes. This result is obvious if we estimate the coupling 
between the torsional motion of a solid and the gravitational field. Using the quadrupole formula, Schumaker 
& Thorne [19] estimate the gravitational- wave damping time scale of the i-modes to be ~ I0 4 years, i.e. 
some fourteen orders of magnitude longer than the typical oscillation time-scale for a 30 Hz mode. The 
upshot of this is that the complex angular frequency w in the standard ansatz for the time-dependence e luJt 
has an imaginary part which is fourteen orders of magnitude smaller than the real part. This is a problem 
for any numerical scheme that aims at determining the complex to. The, perhaps naive, numerical scheme 
we use here is certainly unable to compute the imaginary part of these very slowly damped modes. We 
are, however, able to compute the real part, and hence the frequencies, of the i-modes accurately. The 
obtained results confirm the expectation that these modes can be accurately determined within the Cowling 
approximation [BJ. This is true regardless of whether the superfiuid neutrons are taken into account or not. 
The free neutrons on the other hand do affect the axial mode frequencies, in agreement with the qualitative 
results in |17j . The effect on the global mode frequencies turns out to be at the 10% level, but depends 
strongly on the precise properties of the inner crust. Overall, our results imply that, as far as seismology 
efforts are concerned [BJ, the Cowling approximation should be sufficient. We will return to that problem, 
and detailed implications of the presence of the free neutrons, elsewhere. 

The ii>-modes pose no real technical problems for our numerics and we can determine their frequencies 
and damping times to about 10 digit precision. The results indicate that the inclusion of a more accurate 
treatment of the matter content in the star does not influence the u>-mode spectrum beyond the ~ I0~ 8 
level. This could have been anticipated since the w-modes depend almost entirely on the gravitational field 
of the background star. The background model is hardly affected at all by a more refined matter modelling, 
e.g. in the crust region. 

The plan of the paper is as follows. In section [5] we briefly review the formalism needed for a general 
relativistic treatment of the dynamics of solids coexisting with a (super)fluid. This is followed, in section [3] 
by an explanation of how the free neutrons affect the equations of motion in the linear perturbation regime. 
Our numerical results are then presented in section [U We conclude the paper with section [5] and a discussion 
of the results and possible future developments. The numerical machinery used to solve the equations is 
summarised in an Appendix. 
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2 Formalism 



We begin our analysis by introducing the basic framework used to describe the dynamics of a superfluid 
immersed in a solid lattice. The discussion is based on the theory developed by Carter & Samuelsson 
(CS) [20 , to which we refer for details. We will, however, adapt our discussion to make maximal use of 
the formalism developed in a series of papers by Karlovini et al. dealing with pure solids [211 1221 ESI [24] , 
hereafter referred to as Papers I-IV. We will also, as much as possible, adapt the notation to the recent 
review of relativistic fluids by Andersson & Comer [2S], hereafter AC. 

The theory is derived from a Lagrangian "master function" A which plays the dual role of providing the 
matter part of the total action and also giving the equation of state (EoS). For the purposes of the present 
study it is convenient to decompose the Lagrangian as (see CS) 

A = Ai iq + A sol (1) 

where An q is a function of the particle fluxes n° and n c { , while the clastic properties of the solid lattice are 
contained within A so i. The functional Ai; q is the contribution due to the liquid properties and corresponds 
to the full Lagrangian for a two-fluid system (see e.g. AC). We use the constituent indices 'c' and T to 
distinguish the "confined" baryons in the lattice and the "free" neutrons, respectively. We assume a metric 
signature of the form (— , +, +, +) and use early Latin letters, V, 'd', ... to denote abstract spacetime 
indices, see e.g. [26], omitting 'c' in order to avoid confusion. We will employ liberal index positioning to 
avoid unnecessary cluttering of the formulae. For instance, the total free neutron number density is given 

by 

n 2 = -nfni = -g ab nfn b { (2) 

and similar for the confined baryon number density n c . 

The liquid contribution can be further split into a term that is independent of the relative velocity and 
a piece describing the effects of the relative current, 

Aii q = -p + A ont (3) 

Here we denote the velocity independent term by p since it corresponds to the comoving unsheared energy 
density (which is uniquely defined only for zero relative velocity) as described, e.g. in Paper I. This quantity 
is assumed to describe the minimum energy for a given total baryon density n = n c + rif and we will therefore 
consider it to be a function of n only. The remaining term arises because the state of matter may explicitly 
depend on the relative velocity between the different species of particles in a multi-fluid system. This leads 
to an effect known as entrainmcnt. If only low relative velocities are considered (as will be the case here) it 
is sufficient to consider a slow motion approximation of the entrainment contribution. Since p represents a 
minimum energy state for a given baryon density n it is clear that the leading order term in the entrainment 
term must be quadratic. The question is in what? There are at least two choices. In CS the expansion was 
carried out in terms of the relative flux, , which, geometrically, can be defined locally as the projection of 
the neutron current onto the space orthogonal to the four velocity u a of the solid (or vice versa), 

n\ = h a b n\ , h ab =g ab + u a u b (4) 

On the other hand, in micro-physical calculations it seems customary [571 HE] to use the relative velocity v a 
(or some proxy of this) which is related to the neutron current through 

n? = inf (u a + v a ), u a v a = 0, 1 ={l-v 2 )~ 1 ' 2 , v 2 = v a v a (5) 

Combining and (0 we see that 

v 2 

n\= ins v a => n i = n /T— ^ ^~n 2 f v 2 (l + v 2 ) (6) 

This analysis shows that the two prescriptions agree to 0{v 2 ). However, we also find that the equations 
of motion, which will depend on derivatives of the master function, may differ at 0(v 2 ) depending on the 
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choice of expansion and the level of truncation. In order to be consistent, one should therefore be careful 
to employ the same scheme as the microphysics calculation on which one bases the continuum model. This 
level of consistency is difficult to reach given our present (rather basic) level of understanding, but it sets 
the standard that future work should aspire to. 

In this study we will consider axial linear perturbations of a static background model. The incompressible 
nature of this kind of motion together with the vanishing of the relative velocity in the background implies 
that we will not, even in principle, need to know the master function beyond 0(v 2 ). Thus, for the present 
purposes the problem mentioned above is a non-issue. For more general backgrounds and compressible 
motion (e.g. associated with acoustic oscillation modes) caution is advised. 

Since the microphysics calculations we take our parameters from are performed as an expansion in v this 
is what we will use. Then, making the standard assumption that the liquid is intrinsically isotropic, we can 
write H3 HE] 

Acnt = — m f c (nl f - n c n { ) = n[m { c (^ - 1) = -n f m c w 2 + 0{v 4 ) (7) 

n c 2 

where m c = m f * — m, m f * is the effective (dynamical) mass of the neutrons, m is the neutron mass (which 
we may take to be equal to the proton mass) and n 2 f = —n^n f a . 

The solid contribution is most easily described by an isotropic quasi- Hookean [29, 20] prescription. Then 
it is simply assumed that A so i = —fis 2 where fi = fi{n c ) is the shear modulus and s is a scalar measure of the 
state of strain. Although this description is likely to be adequate in any realistic situation we will nevertheless 
use the more elaborate description discussed in Paper I and only later, in the applications, specialise to 
isotropic quasi-Hookean solids. This will make the resulting equations formally valid for anisotropic solids 
without assuming the Hookean approximation. We will still use a minimal coupling ansatz in the sense that 
we (quite reasonably) assume that the solid contribution is independent of the free neutrons. By using the 
more general matter description, we have the advantage that we can adopt the results of Paper IV more or 
less directly. 

When describing solids it is important to keep track of the reference state relative to which the strain 
is measured. For many simple solids it is possible to do this via a positive definite metric tensor field, the 
matter space metric k a b (see Paper I). One may intuitively think of this tensor as encoding the (3-)geometry 
of the solid (as seen by the solid itself) in an unstrained state. The strain tensor is then given by 

Sab = 7^(hab - n ( T 2/3 fc ab ) (8) 
where h a b is defined in It is advantageous to work with an orthonormal eigenbasis e° of k a b in which 

3 

K b = Y. n VA ( 9 ) 

where we use Greek indices (with no implied sums) to enumerate the basis and are, loosely speaking, 
linear particle densities. In principle we could simply take the solidity contribution to the master function 
to be a function of the n M 's. For practical purposes, however, it is convenient to instead work with the 
eigenvalues a u = -ttj of 

3 

2 /3 

Since the determinant of the matter space metric is given by det(k a b) = n c it is clear that 

n c = nl =1 n^ (11) 

It follows that rfb has a unit determinant so that this tensor only has two independent components. These 
components hold the key information of the material's response to non-compressional distortions and there- 
fore encode the difference from a liquid. Within the eigenvalue formulation it is natural to prescribe the 
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solid contribution as a function of the number density of confined baryons, n c and the parameters a^, 

Asoi = A sol (n c ,a M ) (12) 



making manifest the "minimal coupling" ansatz, i.e. that we assume that the solid's response to deformations 
is independent of the number density of free neutrons. In a quasi-Hookean approximation the solid's function 
of state may be separated in the form 

A so i = fj.(n c )s 2 (a fi ) (13) 

where jj is the shear modulus and s 2 is a scalar measure built from invariants of the strain tensor (see Paper 
I for a discussion). 

To summarise, the above prescription leads to a Lagrangian that can be written 

A = A liq «,<) + A sol (n c ,aJ (14) 
where the standard isotropic two-fluid Lagrangian takes the form 

Anq = -p+ — (ml - m)(nl t - n c m) (15) 

Practically speaking, our matter model is described by three functions of state; the minimum energy density 
p, the effective neutron mass and the solid's function of state which, in the subsequent applications, will 
be encoded in the shear modulus fi. 



2.1 Equations of motion 

The variational procedure described in CS, applied to the Lagrangian density discussed above, leads to a 
stress-energy tensor of the form 

T% = (A - - n d cf i c d )5 a b + n f >6 + « + *\ (16) 

Here 

f _9A_aAuq c _dA_ 

dnf dn a f ' ^ a ~8n- ( ' 

are the momenta of the constituents and ir a b is the (trace-less) solid contribution as derived in Paper I. The 
stress energy tensor may be used as source in Einstein's equations. In addition, we have two equations of 
motion of the form 

2n f a V [a ^, - (18) 
2<V [a ^] + V°7r ab = (19) 

where square brackets indicate anti-symmetrisation. In the following, we will use the fully variational 
description (CS). This means that the constituents are individually conserved, i.e. V a n° = V a rtf = 0. 
The equations (fT8"|) and (fT§|) (corresponding to the Euler equations) together imply the conservation of 
energy momentum V a T a /, = 0. This is, of course, also implied by the Einstein equations via the contracted 
Bianchi identities. Thus the various equations are not independent. This should be familiar from single fluid 
problem where it is well-known that one can opt to work with the Einstein equations alone. For multifluid 
problems more information is required. In the present context it is sufficient to consider a combination of 
the Einstein equations and one of the Euler equations. The simplest choice, since we only need to consider 
the complicated elastic terms once, is to let (|18p do the job. As we will see, the incompressible nature of the 
axial modes together with the fact that the free neutrons are only forced via the entrainment then allows 
us to find the linearised solutions to the mode problem explicitly. We may also comment that the Cowling 
approximation, wherein the gravitational degrees of freedom are neglected, retain the property [6] that it is 
equivalent to consider the weak coupling limit of Einstein's equations together with (fT5)) . Alternatively, one 
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can decide to work only with (fTB)) and (fH?)) . The latter strategy would lead to a problem which is trivially 
related to the corresponding one in Newtonian theory. 

In order to obtain explicit equations we need to determine the momenta. Following the notation of AC 
we write formally 

(20) 
(21) 
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where 

dA 
' Bn\ 

A« := = --^ + 0(, 2 ) « -fmj (24) 

Here we make use of the assumption that p can be treated as a function of n only and denote its derivative 
with respect to n by a prime. The approximate expressions are valid up to 0(v) which is all we will need 
subsequently. We may note here that the generalised pressure (see, e.g. CS or AC) is 

BA BA 
* := A - nX - </4 =np'-p + A sol - n c ^L + {v 2 ) + A sol - n c -^± (25) 

where p denotes the usual (unshcarcd) pressure. Comparing to the result in Paper I we see that \& corresponds 
to the total isotropic pressure defined there. In an unstrained state we have 

-/ P + P _ P _ ^oR^ 
p = w — w m (26) 

n n 

The approximations are valid in the relatively tenuous neutron star crust where p <C p and p « ran. This 
simplification, which is accurate to within ~ 0.1 %, is obviously very useful since we replace a function with 
a known constant. Nevertheless, in our numerical code we use the full expression (which actually means that 
we do not need to assume that the neutrons and protons have the same mass). 

In order to compare the present model to the purely elastic case considered in Paper I it is useful to 
express the energy momentum tensor in a frame adapted to the solid. Formally we then have 

Tab = pu a u b + 2u (a Q b ) + P ab (27) 

where p, Q a and P a b are, respectively, the total energy density, the momentum flow and the pressure tensor, 
all measured in a frame described by the solids four- velocity u a . Putting the pieces together we find that 

p = Pl + 0{v 2 ) (28) 
Q a = Xi ( P + p)v a + 0(v 2 ) (29) 
Pab = ^K b 4 Kab 4 0(v 2 ) = P Q 7 b 4 0(v 2 ) (30) 

where Xf = rif jn is the free neutron fraction and we use the label / to denote the corresponding quantity in 
Paper I. Thus, to order v the free neutrons lead to the presence of a non-zero momentum flow Q a . Obviously, 
a static background is unaffected by this. 



3 Perturbations 

Let us now consider the equations governing axial perturbations around a static background. We base our 
derivation of these equations on the general framework developed by Karlovini }30| and make heavy use of 
the results in Paper IV to which we refer for details. 



G 



The starting point of the analysis is a decomposition of the metric in the form 



g ab =± ab +F~ 1 r ]a r lb , F = rfr) a = (rsin0) 2 , rf ± ab = (31) 

where 77° is the axial Killing vector. Karlovini [3D] showed that the axial part of the metric perturbations 
(denoted here by -f ab ) can be entirely expressed in terms of Srj a so that SF and <5_L a 6 can both be set to 
zero^. The perturbed Einstein equations then take the simple form 

V b {FQ ab ) = K J a (32) 
V a J a = (33) 

where, loosely speaking, Q ab = 2V ' [ a F~ l 8rj b T [ represent the geometric perturbations and 

J a = 28(± ab rj c T bc ) (34) 

encode the matter motion. Hence, it is natural to refer to J a as the matter current. These equations were 
discussed in detail for the case of a single solid in Paper IV. Here we focus on the changes needed to extend 
those equations to the case where the solid is coexisting with a superfluid. As noted above, the superfluid 
degrees of freedom manifest themselves by i) the additional momentum flow Q a in the energy momentum 
tensor, and ii) the need to consider the equation of motion (|18p . In order to work out the explicit perturbation 
equations we found it useful to make use of the eigenvector formulation introduced in Paper I where the 
principal directions of the solid provide an orthonormal tetrad {u a , e°}, where u a is the four- velocity of the 
solid, fi € {1, 2, 3} and we use Greek indices to enumerate the spatial basis vectors. The general expressions 
for the perturbed tetrad were given in Paper IV in terms of the perturbed metrics "f a b — 8g ab and Sk ab which 
in the axial case considered here can be written; 

j ab = 2F -1 77 (a <5?7 b ), Sk ab = 2nlri {a V b )54> (35) 

The quantity cj) is the (inverse) mapping of the azimuthal coordinate on matter space (see Paper IV) and 
we use parentheses to denote symmetrisation. In general, the inclusion of the free neutrons adds four scalar 
degrees of freedom. We take these to be represented by the free neutron density m and the relative velocity v a 
(which, due to the constraint u a v a = has three degrees of freedom) . For the axial case it is straightforward 
to show that we must have (for a static background) 

<5rt c = Sri{ = 8n c f = (36) 

It follows that 

8B C = 5B l = SA ci = (37) 
so that the perturbed momenta are given by 

8nl = B c 6n c a + A cf 8nl (38) 
S^^B'Sni + A^Snl (39) 

For the perturbations we shall consider the relative velocity to be aligned with the axial Killing vector. 
Therefore we may write it in the form 

Sv a = v a = vrf (40) 

where we have dropped the perturbation symbol S since no confusion can arise. This leads to a simple 
expression for the perturbed neutron momentum, 

5^ a = [{p' + m{)v-p'u b K b ]i la (41) 



1 In order to avoid introducing yet another fj, we have chosen to express the metric perturbations in terms of rj a rather than 
[i a = F~ 1 r] a used in |30| and Paper IV. Note that &r] a = F6fj, a 7^ g a b^V a = where the latter equality is due to a partial gauge 
fixing (which is just for convenience since the present formalism is gauge invariant). 
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where K a = V ' a 5<j> — F~ 1 5rj a is a gauge invariant one-form. As discussed in paper IV, it encodes the complete 
nature of the perturbations of the solid. On a static background, and only considering the non-stationary 
(oscillatory) perturbations of equation (fT5|) we find the remarkably simple result 

SfA = (42) 



This implies, via (|41|) , that we have an algebraic solution for the perturbed relative velocity 



u a K a ^—u a K a (43) 
p' + ml m { 

Note also that this shows that the vorticity remains zero to the first order accuracy considered here. Thus, 
no vortices are created or destroyed dynamically to this order. This simple result may seem surprising at 
first, but it is a direct consequence of the incompressible nature of axial perturbations. The variational 
principle we use is based on varying the flowlines of the particles. Incompressibility implies that the total 
number of flowlines in a 3-volume is conserved and hence restricts the variations further to the extent that 
only a single displacement vector will be needed to describe the motion of the full system. 

Armed with the solution to the equation of motion for the free neutrons we are in a position to evaluate 
the effect the free neutrons have on the perturbed energy momentum tensor. As already pointed out, the 
only difference from the purely elastic case is the presence of the momentum flow Q a , Perturbing ([%9"|) and 
inserting the result in (|3"4"1) we readily find 

J a = 2(p + Pt )FS ab K b (44) 

where p and p t are the total energy density and the tangential pressure, respectively, and the difference 
induced by the superfluid is the modification of the shear wave "metric" : 

S ab = S ab + xi — ~ r £±£u u 6 (45) 
p' + m{ P + Pt 

Here 

cab a b , 2 a b , 2 a b Iag\ 

b = —u u +v r± e 1 e% +v tJ _ t e 2 e 2 (46) 

is the corresponding metric for purely elastic matter expressed in terms of the eigenvector basis e^ a of the 
solid lattice. Meanwhile, v r ± and v t ±t are the shear wave velocities in the radial and tangential directions, 
respectively (see Paper I). Hence, the effect of the neutrons is entirely contained in the factor 

X = l-*f / f ' P+ J (47) 
p' + mi P + Pt 

From Paper I we know that the change of the pressure and density due to anisotropy is of the order ~ (is 2 . 
In the bulk of the crust we have the ordering p, <C p <C p and in addition s 2 < 0.1 [31J . Hence we expect that 
the approximations p ~ p and pt ~ p will always hold in realistic neutron star crusts. They will, of course, 
be exactly true on an unstrained background. We therefore note that since p' ~ m we have (to ~ 0.1% 
precision) 

Tfl 

X«l-x f — (48) 
m\ 

in perfect agreement with result from a recent analysis of incompressible plane waves in the corresponding 
Newtonian problem [17j . Thus, to "upgrade" the equations in Paper IV we only need to insert the factor x 
appropriately. The general equations for a static background and I > 2 thus become 

-W t +W r - ? -W r - e 2l/ 4^ = (49) 

-W r + W t + -Wt + e 2v ^±!-^ = Q (50) 

-Lrip + Ev 2 ^ 2 ^- 1 ^) 1 + (£v 2 ± + L)rW t = (51) 
-X£r<p + L(rip)' - (x£ + L)rW r = (52) 
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cf. equations (98)-(101) in Paper IV. Here r is the radial measure in Schwarzschild coordinates in which the 
line element takes the form 



ds 2 = -e^di 
We are also using L = (I — + 2), 



IVA-,1 + e 2A d? ,2 + r 2( d6 ,2 + sin 2 

£ = 2 K r 2 {p + Pt ) 



(53) 
(54) 



and dots and primes refer to derivatives with respect to time and the Regge- Wheeler radial coordinate r* 
which is given by 



dr* = 



(55) 



The total energy density p and the tangential pressure pt (see Paper I) have been retained to keep the formulae 
valid for anisotropic backgrounds. We will later set them to the isotropic values p and p, respectively. For 
an explanation of the dependent variables, see Paper IV. We see that only the last equation is changed by 
the presence of the free neutrons. 

For completeness we also give the modified wave equations (of which only the second, describing the 
shear waves, is modified) 



-W t + WJ + e 



6m 1 



£v? ± 



1(1 + 1] 



-X ( P + 



£r 
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[£(x-v^)rmy x£' 

£r £ 



(56) 
(57) 



The boundary conditions, namely that Wj, W r and tp are everywhere continuous, do not change. Thus, we 
now have all the information needed to solve the axial oscillation problem including the free neutrons. 



4 Results 

The basic strategy for solving the perturbation equations ([4*9")) - ("5"2")) . in the case where the free neutrons 
are ignored, (\ — 1) has already been outlined in Paper IV. Given that the free neutrons do not change the 
problem formally, the same strategy can be used in the more general setting discussed here. Nevertheless, 
it is worthwhile providing some detail on the methods we have used to to solve the equations. Interested 
readers can find a discussion of the relevant points in the Appendix. 

We build our background neutron star models using realistic tabulated equations of state, see table [T] 
The two older equation of state tables "A" [32] and "B" [33] were taken from the distribution of the rns 
code, see e.g. [3(5], whereas "APR" [31] was provided by G. Ravenhall. In the crust the equation of state is 
taken from Douchin & Haensel [37] (DH jf] and we use the shear modulus for a Coulomb lattice as calculated 
by Ogata & Ichimaru [39jf|. The effective mass needed for the entrainment was taken from Chamel [5T] . We 
emphasise that, since his calculations were performed for a different EoS than the ones we use, our model 
for the effective mass is inconsistent. In addition, we are only aware of data for four particular densities. 
Since we need data for all crust densities, we use an analytic expression that approximately passes through 
the available data points, see [17]. This is, obviously, not satisfactory but it is the best that we can do at 
the present time. Further work on the properties of the crust "beyond" the EoS (minimum energy state) 
should be encouraged. Even though our model is somewhat ad hoc, it should provide insights into the basic 
dynamics of the problem. 

As a first test of the numerical code we determined the w-modes for a number of cases: 

2 In fact, in the outer crust we employ the results described in |38| using more recent measurements for the binding energy 
of nuclei, but the end result is practically indistinguishable from DH. 

3 The effects of using the very recent results of Horowitz & Hughto [40] . obtained from molecular dynamics simulations, will 
be examined in a forthcoming paper. 
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Model 


EoS 


Po (10 15 g/cm 3 ) 


M(M P) ) 


R (km) 


M cor JMp,) 


Rrr,rr fklll) 

"LUI (J \ / 


8 


Al(a) 


A 


1.259 


1.04959 


9.89188 


- 


- 


0.156751 


Al(b) 


A + DH 


1.259 


1.04806 


9.90598 


1.03690 


8.99246 


0.156300 


A2(a) 


A 


4.110 


1.65241 


8.37048 


- 


- 


0.291632 


A2(b) 


A + DH 


4.110 


1.65019 


8.37450 


1.64746 


8.10619 


0.291100 


Bl(a) 


B 


1.995 


0.97068 


8.76793 


- 


- 


0.163549 


Bl(b) 


B + DH 


1.995 


0.97045 


8.78057 


0.96302 


8.01716 


0.163274 


B2(a) 


B 


5.910 


1.41131 


7.07165 


- 


- 


0.294828 


B2(b) 


B + DH 


5.910 


1.41003 


7.07679 


1.40843 


6.85589 


0.294348 


APRl(a) 


APR 


0.750 


0.91989 


11.5963 


- 


- 


0.117189 


APRl(b) 


APR + DH 


0.750 


0.92026 


11.7203 


0.89627 


10.1619 


0.115995 


APR2(a) 


APR 


2.750 


2.19428 


10.0059 






0.323969 


APR2(b) 


APR + DH 


2.750 


2.19436 


10.0238 


2.19071 


9.77753 


0.323404 



Table 1: Background models used in our study. We consider three equations of state, labelled A [32], B [33J 
and APR 34J. The reason for this particular choice is simply that they were previously studied by BBF 
[35] . This allows us to compare the results for the w-modes directly. For each equation of state we consider 
the same two central densities as BBF. For each of these models we also substitute the equation of state by 
DH in the crust. As the data in the table shows, this leads to slightly different total masses and radii. For 
each model we provide the central density po, total mass M, total radius R and compactness 8 = M/R. For 
the models where the crust-core transition pressure is known (i.e. for the DH crust) we also give the mass 
M COTe and radius R corc of the core. 

(la) using the fluid equations with the available tables, as described above, all the way to the surface, 
(lb) as (la) but substituting using the EoS of Douchin & Haensel in the crust region 

(2) artificially setting the shear modulus to zero, 

(3) using the full equations but artificially setting x = 1, and finally, 

(4) using the full equations including the model for the superfluid neutrons. 

Case (la) can be directly compared to the results of Benhar et al. 35j (BBF). We compare a selection 
of typical results in table (2] This comparison shows that the results are generally in good agreement. The 
small differences could be due to slightly different tabulations of the EoS or different interpolation schemes 
(we use logarithmic interpolation). We also expect our treatment of the surface of the background model 
to be more accurate. To rule out the possibility that the difference is due to a bug in our code we also 
calculated w-modes for a polytropic equation of state and compared to the results of Andersson & Kokkotas 
[32]. In this case (where the ambiguity in the treatment of the EoS is eliminated) we find perfect agreement. 
A third test of the code is provided by comparing cases (lb) and (2). The agreement is excellent also in this 
case. 

In order to investigate what effect the elastic solid and the superfluid neutrons have on the w-modes it 
is relevant to compare cases (2)-(4). Since the w-modes are of gravitational origin and the matter motion 
is limited to the relatively tenuous crust we expect a very small effect. The numerical results confirm this 
expectation. We find that the relative influence of the crust and the superfluid is less than ~ 10~ 8 (see table 
[3] for some typical results). This difference is roughly at the same level as the expected accuracy of our code. 
Thus, for all practical purposes, the axial w-modes are unaffected by the crust properties. 

Turning to the torsional i-modes, where one would expect the superfluid to have a significant effect, we 
compare cases (2) and (3) to our previous results [S], that were obtained within the Cowling approximation. 
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Model / (a) (kHz) f {b) (kHz) / B bf (kHz) r (a) {(is) r {b) {(is) t B bf {(Js) 



Al 


q 701 1 oe 




Q 7fi 


91 574378 


21 540601 


21.6 


A 2 


q 1 07043 


Q 1 1 91 39 


Q 1 1 


70 028558 


71 577816 


72 4 


Bl 


11.24895 


11.25115 


11.2 


20.003988 


19.999529 


20.2 


B2 


10.65463 


10.66564 


10.6 


70.601159 


70.188532 


71.7 


APR1 


8.865023 


8.865084 


8.82 


19.195766 


19.203608 


19.4 


APR2 


6.725440 


6.725176 


6.69 


158.52087 


158.53555 


165.3 



Table 2: The data in this table compares the w-modes that we have determined to the results of BBF. 
We give the frequencies and damping times for the lowest I = 2 iu-modes for completely fluid stars. The 
subscripts refer to the cases discussed in the text, i.e. (a) uses the tables discussed in the main text and 
the fluid equations, (6) substitutes the DH EoS in the crust and use the solid equations with (i set to zero. 
Finally, BBF refers to the results in [55] . 



Model 


/ (kHz) 


r (ms) 


Fluid 


9.11913825 


0.0715778249 


Solid 


9.11913834 


0.0715778266 


Solid + SF 


9.11913833 


0.0715778264 



Table 3: The data in this table illustrates the effects of the matter model. We give the frequency and damping 
time of the lowest I = 2 w-mode for the background model A2(b) (see table [T]) for different matter models. 
The three models are: a completely fluid star {fi — 0, 'Fluid'), a star with a solid crust, but neglecting the 
free neutrons ('Solid') and the complete model including the neutrons ('Solid +SF'). The results show that 
the effects due to the more detailed treatment of the crust physics are tiny. The results have converged (after 
increasing the resolution) to the stated precision, but one should be aware of the fact that the roundoff error 
in the calculation is, as discussed in the main text, of the order of ~ 10~ 10 . 

Since the estimated (via the quadrupole formula) gravitational-wave damping tim^] for the i-modes is about 
t ~ 10 4 years [H] we expect the real part of the frequency to be well approximated by the Cowling results 
and the imaginary part to be very small. Indeed, a typical I = 2 fundamental torsional mode has / ~ 30 Hz 
so that 

$t(w) = 2tt/ - 1.9 x 10 2 s" 1 (58) 

In contrast, the imaginary part is 

3f(w) = - ~ 0.8 x 10~ 12 s" 1 (59) 

T 

some fourteen orders of magnitude smaller than the real part. Since our numerical code uses an adaptive 
step Runge-Kutta integrator we can control the accuracy in each step. However, using a small error tolerance 
per step induces an uncontrollable truncation error (we use double precision so that the truncation error is 
~ 10~ 15 ) which, if turning down the tolerance on the error in each step, requires more steps to complete the 
integration. A simple analysis, using the de facto used steps together with the prescribed tolerance, suggests 
that our code cannot reach an accuracy beyond about 1 part in 10 10 . Thus, we should not expect to be 
able to determine the the damping time of the t-modes directly. The numerical calculations confirm this 

4 Note that the estimates in |19] concern the damping time for the energy, while we need the damping time associated with 
the amplitude. However, given the order of magnitude nature of our discussion of this point we do not bother with the different 
factors of 2 that relate these quantities. 



11 



Model 




2/0 


(Hz) 


2/1 


(Hz) 


2/2 


(Hz) 


Alibi 


31 


097 


(31.069) 


881.15 


(881.15) 


1449.2 


(1449.4) 


Al(b) + SF 


33 


897 


(33.872) 


957.93 


(957.92) 


1537.7 


(1537.8) 


A2(b) 


27 


169 


(27.160) 


1794.2 


(1794.2) 


2962.9 


(2963.3) 


A2(b) + SF 


29 


672 


(29.665) 


1950.3 


(1950.2) 


3164.7 


(3165.1) 




34 


549 


(34.525) 


1031.5 


(1031.5) 


1697.3 


(1697.5) 


Bl(b) + SF 


37 


666 


(37.645) 


1121.4 


(1121.4) 


1802.1 


(1802.2) 


B2(b) 


31 


874 


(31.867) 


2144.8 


(2144.8) 


3542.1 


(3542.5) 


B2(b) + SF 


34 


812 


(34.806) 


2331.4 


(2331.4) 


3783.6 


(3784.0) 


APRl(b) 


28 


887 


(28.841) 


584.00 


(584.00) 


955.60 


(955.66) 


APRl(b) + SF 


31 


449 


(31.408) 


634.90 


(634.89) 


1008.0 


(1008.1) 


APR2(b) 


20 


739 


(20.731) 


1649.4 


(1649.4) 


2724.7 


(2725.0) 


APR2(b) + SF 


22 


655 


(22.647) 


1792.9 


(1792.9) 


2912.1 


(2912.5) 



Table 4: This table provides a sample of results for the first few quadrupole (/ = 2) t-modes. The background 
models are explained in table [T] and the extra label "+ SF" indicate that we have taken the free neutrons 
into account. The results using the Cowling approximation are shown within parenthesis for comparison. 
Since the code we use for the Cowling approximation is optimised for speed rather than accuracy it can 
only deliver 5-6 digit precision which is why we only quote this accuracy here. The calculations reported in 
this work is much higher, but this is most likely irrelevant for astrophysical applications. We note that the 
Cowling approximation in the cases tested is accurate to better than about 0.01 %. We also note that the 
case where the free neutrons are taken into account typically gives ~ 10 % higher frequencies, but we stress 
that this is model dependent and the effects can be much larger (or smaller). 

expectation. When we feed the output of the integrator into the Miiller root-solver we find that we cannot 
determine the imaginary part of the mode frequency. The real part, on the other hand, converges nicely 
to more than ten digit precision (while the imaginary part oscillates around zero). For this reason we are 
confident that we determine the oscillation frequency accurately. Typical results are provided in tabled In 
the table, we compare the results obtained using the Cowling approximation (the case where the superfluid 
is neglected was discussed in [5] and the analysis of the problem including the superfluid is currently being 
finalised). 

5 Discussion 

In this paper we have considered the global axial quasi-normal modes of relativistic stars with a (possibly) 
superfluid core and a solid crust penetrated by a superfluid component. Our results provide the first detailed 
analysis of this problem in general relativity, and represent a key step towards the modelling of realistic 
neutron star dynamics. In fact, apart from the linear approximation common to all mode studies and 
the omission of the magnetic field, our treatment does not impose any significant approximation^. The 
discussion does, however, highlight the need for improved microphysics models. While we have made an 
effort to use models that are as "realistic" as possible, it is clear that the input parameters that we have used 
are somewhat inconsistent. This is entirely due to the lack of complete data from microphysics studies. Future 

5 The treatment of the solid component assumes a conformally deforming solid, see Paper I, but as this class incorporates, 
e.g. cubic symmetric lattices and isotropic solids we do not think that this constraint imposes any important restrictions at the 
present time. 



12 



tabulated equations of state need to provide, in particular, the superfluid entrainment parameters (both in 
the crust and in the core). Once such data becomes available, it will be straightforward to incorporate it in 
our computational framework. The issue concerning the neutron star magnetic field is more challenging. This 
problem will require serious thought, especially if we want to account for the presence of superconducting 
protons in the core. 

Specialising to an isotropic solid we determined both the axial gravitational w-modes and the elastic 
torsional i-modes. In the case of the w-modes we calculated both the frequency and the damping time with 
high precision. The obtained results confirm the expectation that these modes are very weakly influenced 
by the presence of a solid/superfluid component. We also confirmed that the frequencies of the t-modes are 
only weakly affected by the coupling to the gravitational degrees of freedom. We were, however, unable to 
directly determine the damping time for these modes. This is likely to be completely irrelevant from an 
astrophysical point of view, since other mechanisms will dominate the damping of these modes. Of course, 
from an academic point of view the situation is not entirely satisfactory. As a matter of principle, one would 
like to have a direct determination of the damping times. Of course, a direct application of the quadrupolc 
formula should give reliable results. A possible future option would be to take advantage of the enormous 
difference between the real and imaginary parts of the quantities that appear in the perturbation equations 
and perform a "second perturbation" . One can easily check that this leads to a set of decoupled equations 
for the real parts and a set of equations for the imaginary parts sourced by the real parts. We have not yet 
tried to implement such a scheme. 

As far as any magnetar seismology analysis is concerned [6] , the main lesson to learn from our analysis is 
that one does not need an accurate treatment of the gravitational degrees of freedom. However, the presence 
of the superfluid component in the crust is important. We are currently investigating this problem in more 
detail within the relativistic Cowling approximation, and hope to be able to report on the results in the near 
future. Having developed the framework for the axial perturbations for relativistic neutron star models with 
elastic and superfluid components, we also need to consider the (generally more complex) problem of polar 
perturbations. Developments in this direction are also under way. 
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A Appendix: Numerical formalism 

In this appendix we discuss the methods we used to solve the axial perturbation equations - (f5^|) . In 
paper IV the problem was analysed for the case x — lj *- e - when the free neutrons are ignored. Given that 
there is no real formal difference between the two systems that need to be solved, we can use the same 
strategy also for this more general problem. 

In the interior of the star we use a standard Runge-Kutta ODE solver from the GSL libraries [43J. The 
exterior perturbations are determined by a pseudo-spectral method (see [H] and below). 

A.l Background solution 

A static, spherically symmetric, unstrained solution does not depend on the elastic/superfluid nature of the 
matter content since both the strain and entrainment contributions vanish. Thus, the background solution 
is, in this case, identical to that of a perfect fluid. It is well known that the use of a radial coordinate as 
independent variable, as in the standard TOV equations, leads to numerical difficulties near the surface. 
One problem is the steep gradient of the fluid variables near the surface and another is that the surface itself 
is determined by the condition that the radial pressure vanishes. That is, it is determined by one of the 



13 



dependent variables. Lindblom [35] suggested that a better strategy, as far as the background is concerned, 
is to define the two dependent variables 



u = r , v = — (60) 
r 



and use the relativistic enthalpy h defined by 



dh = (61) 
P + P 



as the independent variable. Then the equations of structure take the form 

du 4u(l - 2v) 



dh Kpu + 2v 



: u' (62) 



dw , „ x Kpu — 2v u' „ . ,„ . 

— = -{l-2v) — — = —Upu-2v) (63 



In order to ensure a regular centre the integration is started at finite radius using the expansion 

u= 12(h -h) 
k(pq + 3p ) 
2p (h -h) 

V = h • • • 

Po + 3p 

where a subscript '0' will be used throughout our discussion to denote evaluation at the centre. 



(64) 
(65) 



A. 2 Perturbations 

As we are interested in quasi-normal modes we make the standard assumption that the time dependence 
is given by exp(iwt) where w is a complex constant whose real and imaginary parts represent the angular 
frequency and the damping timescale, respectively. Inside the star we solve the perturbation equations 
together with the background equations. 



A. 2.1 The vacuum region 

In the surrounding vacuum region the axial perturbation equations reduce to the standard Regge- Wheeler 
equation. For quasi-normal modes we require that the solutions are outgoing waves at null infinity. This 
problem was discussed in detail by Samuelsson, Andersson and Maniopoulou [2] and we apply their code as 
it is. Then, given the mass M, radius R and angular frequency to we obtain the surface value g s of a certain 
phase function 

" U (66) 



9\r=R 



1 dip 
ip dr 



-R 



A. 2. 2 The fluid core 

Following Paper IV we adapt the independent variables to the centre of the star by defining 

Xi = r-^Wt (67) 

X 2 = -iu)e h °- v °r- l W r (68) 

X 3 = -ioje ha - h r- l - 1 iP (69) 

X 4 = u 2 e 2h °- h -"T- l ip (70) 
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where v s = i In (l — denotes the surface value of v. In the fluid region the equations then reduce to 



dA"i _ u' 
~dh ~~ 2i7 
dX 2 _ vf_ 
~dh ~~ ~2ii 
and the constraints 



^h — hn 



(I + 2)Xx + 



-.Xo 



?h-h 



y/T=2v 



(/ + 2)(i- 1) 



_ 1^2(^0-/1) _ ,, .2-2(A -f,), 



Afj + (Z - 1)^ 



(71) 
(72) 



X 3 = -e h °- h X 1 , 



Xi = -e ll °-' l X 2 (73) 

We may note here that the perturbations in the core do not depend on the multifluid nature of the medium 
{i.e. x does not appear). This is due to the fact that, like an ordinary perfect fluid, this type of matter 
cannot sustain axial oscillations, but only stationary currents. Hence, the only oscillations that remain are 
those associated with the gravitational degrees of freedom (the w-modes). This means that, one cannot use 
the axial w-modes to distinguish between non-rotating single- and multi-fluid stars. If the star is rotating, 
however, it is known that the "axial-led" inertial modes will depend on the superfluid nature of matter (46] . 
In order to make sure that we have a regular solution near the origin we expand the solution according 

to 



, 6e 2 ^-^ 2 + K (l + 2)[3p - (21 - l)p } 
1 — — — — — (Hq — n) + 



X 9 



Xi -(1 + 2) 



(2l + 3)n(po + 3p ) 
6(1 + 4)e 2 ( fe »-^a; 2 - (I + 2)(l - l) K [3po + (21 + 7)p ] 
(21 + 3)n(p + 3p ) 



(h -h) + 



(74) 
(75) 



where Xi represent the arbitrary scaling of the solutions and one should keep in mind that we are allowed 
to rescale the solutions at our convenience. We check that the solution does not depend on the (small) value 
chosen for ho — h. 



A.2.3 The crust 

In the crust we have chosen to integrate a slightly different set of equations with dependent variables given 

by 



yi = rW t 

y, 2 = _L Wr 

icvr 

34 = r(W t - M/0 



34 = -<p- 

r 



1 



-W r 



(76) 
(77) 
(78) 
(79) 



For these variables 34-4 are constrained to vanish in a fluid and 34 = also in vacuum. Moreover, 34-3 are 
everywhere continuous. Since we lack concrete information about any anisotropy in the solid we specialise 
the equations to the isotropic case where 



v?x = Vt Xt 



I 1 



P + P 



P = P, 



Pt = Pr = P 



The set of equations we solve then becomes 



yi,h = -p« 
p 

y^,h - — r— r 



y s , h 
y^.h 



2kPu r 



L 



2nu 2 uj 2 jl 



tu 2 y 2 - 2ne 2 ^~ h ^(y 2 - 34) 

\u 2 uy 1 -Le 2 ^- h \y 1 -y 3 ) 
X uj 2 u(p + p)y A + Le 2 ^^p(y 2 - y 4 ) 

u 2 y 3 + 2ne 2 ^' h ^p(y l - 34) 



LP 



(80) 

(81) 
(82) 
(83) 
(84) 
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where P = u' 'e u ^ h /2y/u(l - 2v). 



A. 2.4 Boundary conditions 

As outlined above, the boundary conditions both at the centre of the star and at infinity are already taken 
care of in our scheme. It thus remains to ensure that the proper jump conditions are satisfied at the top 
and bottom of the crust. At these interfaces we must ensure that 3T-3 are continuous. As discussed in 
Paper IV, we have no local conditions on the fourth variable. However, a generic choice of 3^4 at the inner 
(say) boundary will not admit a continuous solution at the outer boundary. For this reason we solve for 
two linearly independent solutions in order to find the unique^ linear combination that satisfies all boundary 
conditions. 

There is one difficulty with this approach, however. The first term in equation (|84p is proportional 
to Near the surface of the star the boundary conditions require that 3^3 — > which balances the 

smallness of the shear modulus jl in that region. Unfortunately, neither of the two outgoing solutions that 
we compute would typically satisfy the boundary conditions at the surface (only a specific linear combination 
will). Hence, this term will typically become very large close to the surface which makes it difficult to find 
an accurate solution. To remedy this situation we integrate the eigenfunctions in the crust both from the 
top and the bottom to some intermediate matching point. We check that the choice of matching point does 
not affect the solutions. 

The explicit matching conditions are as follows. At the fluid-solid interface we have 

yi = r< I+2 >Ai (85) 

y 2 = e v °rV-Vcj- 2 X 2 (86) 
3>3 = (87) 
34 = free 



where the Xi's are given by the fluid solution. We define the two linearly independent solutions with the 
interface values 

jf» = [ r ( |+ ^ 1 ,e^(^ W - 2 Z 2 ,0,0], yf } = [0,0,0,1] (89) 

which guarantees (using the freedom to rescale the fluid solution) that the linear combination 3^ tot ^ = 

Ci^j + C^yf 2 satisfies the boundary conditions at the surface for any values of d. At the surface we 
arbitrarily fix the overall normalisation of the solution by setting y± — 1 where we will use a tilde to 
distinguish the outer crust solution from the inner. From the vacuum value of the phase function (|66p we 
then obtain the boundary value of 3^2 , 

y2 R\P (90) 

In addition we require that 3^3 = 0. Thus, we choose the independent solutions via the interface values 

3? } = [l,3Mff,),0,0], 3>f' = [0,0,0,1] (91) 

This again guarantees that the linear combination y^ ^ = y^ 1 ' + C^yf"' satisfies the boundary conditions 
at the interface for any value of C3 . Note that we only have a single coefficient in the linear combination due 
to the choice of normalisation. 



6 The solution will be unique only if one assumes that ^4 is continuous in the crust. This would seem natural, but is in fact 
not required by the equations. One could imagine discontinuities at e.g. phase transitions between different regions in the solid. 
We take to be continuous, corresponding to the assumption that the crust behaves as a single solid. It is our understanding 
that the sharp discontinuities found in calculations of nuclear matter at absolute zero temperature are likely to be smoothed 
in a real neutron star. At least on the length-scales of interest in a global mode analysis. Moreover we find it highly unlikely 
that, even if such discontinuities appear, the layers will slip freely along the boundaries. 
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We have now ensured that the boundary conditions are fulfilled everywhere except at the matching point. 
Here we require that all 3^ are continuous. Thus we demand that 

Ci^ 1} + c 2 yl 2) = yP + c 3 yl 2) (92) 

This condition can be put in matrix form as 

Mi; = (93) 

where 



M = 













\ 










y? 










yp 






V 


yP 


yP 


yP 


*i 2) 


/ 



and 

w = [Ci,C 2 ,-l,-C 3 ] T (94) 

In order for a solution exist we need det(M) = det[M(w)] = which is the eigenvalue equation that gives 
the quasi-normal frequencies. We solve the mode condition using a standard complex root finder based on 
Miiller's method. 
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